In [1]:
from functools import partial
from pathlib import Path
In [2]:
import pandas as pd
import numpy as np
In [3]:
import plotly.offline as py
py.init_notebook_mode()
In [4]:
import holoviews as hv
hv.extension('plotly')
In [5]:
import matplotlib.pyplot as plt
In [6]:
%matplotlib inline
In [7]:
# first value is default
# rest are alternative value for each realization
PARAMS = {
'STELLAR_BARYON_FRAC': (0.05, 0.001, 1.),
'STELLAR_BARYON_PL': (0.5, -0.5, 1.),
'ESC_FRAC': (0.01, 0.001, 1.),
'ESC_PL': (-0.5, -1., 0.5),
'M_TURNOVER': (5e8, 1e8, 1e10),
't_STAR': (0.5, 0., 1.),
'L_X': (40.5, 38., 42.)
}
In [8]:
from dautil.plot import iplot_column_slider, plot_column_slider
In [9]:
def read_global_evolution(path):
df = pd.read_csv(path, delim_whitespace=True, header=None, names=('zp', 'filling_factor_of_HI_zp', 'Tk_ave', 'x_e_ave', 'Ts_ave', 'T_cmb*(1+zp)', 'J_alpha_ave', 'xalpha_ave', 'Xheat_ave', 'Xion_ave'))
df.set_index('zp', inplace=True)
return df
In [10]:
def parse_case(string):
for level, idx in enumerate(map(int, string.split('_'))):
if idx != 0:
break
if idx == 0:
return 'default'
else:
param = list(PARAMS)[level]
return '{}={:.6}'.format(param, PARAMS[param][idx])
In [11]:
basedir = Path('~/21cmfast/21cmFAST/Parameter_spaces').expanduser()
In [28]:
df_path = pd.DataFrame((basedir.glob('**/global*')), columns=['path'])
In [30]:
df_path['case'] = df_path.path.map(lambda path: path.parent.parent.parent.name)
In [32]:
df_path['name'] = df_path.case.map(parse_case)
In [35]:
# filter out bad cases
df_path = df_path[~df_path.case.isin(('0_0_0_0_0_1_0',))]
In [36]:
df_path
Out[36]:
In [65]:
df = pd.concat(
(read_global_evolution(path).x_e_ave for path in df_path.path),
axis=1,
# keys=df_path.name
keys=df_path.case
)
In [39]:
df.sort_index(inplace=True)
In [ ]:
df.to_hdf(
'21cmfast-x_e_ave.hdf5',
'df',
format='table',
complevel=9,
fletcher32=True
)
In [13]:
py.iplot(iplot_column_slider(df))
In [14]:
py.plot(iplot_column_slider(df), filename='x_e.html')
Out[14]:
In [53]:
df.head()
Out[53]:
would be z values when x_e crosses 1 if interpolated linearly
In [62]:
x_0 = df.index[0]
dx = df.index[1] - x_0
y_0 = df.iloc[0]
dy = df.iloc[1] - y_0
(1. - y_0) * dx / dy + x_0
Out[62]:
Generate class .ini files¶
In [117]:
INI = '''H0 = 67.32117
omega_b = 0.02238280
N_ur = 2.03066666667
omega_cdm = 0.1201075
N_ncdm = 1
omega_ncdm = 0.0006451439
YHe = 0.2454006
tau_reio = 0.05430842
n_s = 0.9660499
A_s = 2.100549e-09
non linear = halofit
output = tCl,pCl,lCl,mPk
lensing = yes
root = output/{case}-
write warnings = yes
write parameters = yes
input_verbose = 1
background_verbose = 1
thermodynamics_verbose = 1
perturbations_verbose = 1
transfer_verbose = 1
primordial_verbose = 1
spectra_verbose = 1
nonlinear_verbose = 1
lensing_verbose = 1
output_verbose = 1
format = camb
reio_parametrization = reio_inter
reio_inter_num = {reio_inter_num}
reio_inter_z = {reio_inter_z}
reio_inter_xe = {reio_inter_xe}
write thermodynamics = yes
'''
Binned¶
In [66]:
def binning(array, binwidth=4):
n, m = array.shape
n_trunc = (n // binwidth) * binwidth
if n_trunc != n:
print(f'array truncated from {n} to {n_trunc}.')
array = array[:n_trunc]
return array.T.reshape(m, -1, binwidth).mean(axis=-1)
In [105]:
def format_float_array_trunc(array):
string = ','.join(f'{i:.15}' for i in array)
print(len(string))
return string
In [83]:
def format_float_array_st(array):
string = ','.join(f'{i:.16E}' for i in array)
print(len(string))
return string
In [81]:
def format_float_array(array):
string = ','.join(map(str, array))
print(len(string))
return string
Unbinned¶
In [67]:
# setting boundary values according to documentation on reio_inter
df.loc[0.] = -2.
df.loc[6.] = -1.
df.loc[30.] = 0.
In [69]:
df.sort_index(inplace=True)
In [71]:
df.head()
Out[71]:
In [82]:
Out[82]:
In [74]:
for name, series in df.items():
break
In [84]:
series.head()
Out[84]:
In [106]:
for name, series in df.items():
assert len(format_float_array_trunc(series.values)) < 1024
In [123]:
outbasedir= Path('~/21cmfast/21cmFAST/Parameter_spaces').expanduser()
In [124]:
outbasedir
Out[124]:
In [125]:
reio_inter_num = df.shape[0]
reio_inter_z = format_float_array(df.index.values)
for case, series in df.items():
reio_inter_xe = format_float_array_trunc(series.values)
with open(outbasedir / f'{case}.ini', 'w') as f:
print(INI.format(case=case, reio_inter_num=reio_inter_num, reio_inter_z=reio_inter_z, reio_inter_xe=reio_inter_xe), file=f)
In [122]:
for case in df:
(outbasedir / f'{case}.ini').unlink()